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【摘要 】 极 端 不 平衡 数据 定义 为 上 


变量 或 因 


变量 指标 的 取 值 
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现 严重 比例 失衡 的 数据 ， 例 


如 病例 -对 照 极 度 不 平衡 、 疾 病 发 病 率 极 低 、 生 存 数 据 大 量 删 失 以 及 遗传 位 点 为 低频 或 罕见 


变异 等 。 在 此 情境 下 ， 
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HER. Cox 比例 风险 模型 等 参数 


偏离 正 态 


以 控制 


类 错误 。 


近年 来 ， 随 着 超大 型 人 群 队列 全 基因 


恨 设 检验 的 经 典 统计 量 


组 关联 研究 资源 


的 日 益 共 
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近年 来 ， 随 着 大 型 生物 村 
科 气 的 机 遇 与 挑战 并 存 。 一 方面 ， 更 加 充足 的 村 


Eo 高 效 准 而 
突出 。 为 此 ， 本 文系 统 地 进行 了 方法 学 概述 。 
阐述 极端 不 平衡 数据 对 统计 
正方 法 : Firth 校正 和 毅 点 近似 方法 ， 最 后 ， 简 介 极 端 不 平衡 基因 
本 文 为 极端 不 平衡 数据 的 统计 分 析 提 供 理 论 参 考 和 应 


因 组 关联 研究 ， 极 前 


fj Ach FEE 


记 物 ， 男 一 方面 ， 采 


导致 严重 的 一 类 错误 


ESA 
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常规 统计 方法 
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并 或 非 独立 样本 极端 不 平衡 数据 的 统计 需求 日 益 
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m2 


分 布 的 影响 ， 然 后 


首先 ， 综 述 
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不 平衡 数据 ，Firth 校正 ; 


推荐 。 


本 库 和 全 基因 组 测序 数据 


分 析 极 端 不 平衡 数 
组 关联 在 


的 日 益 释 放 与 共享 ， 遗 传 数据 深度 


E 
NJ 


本 和 测序 深度 有 助 于 识别 新 的 疾病 相关 标 


(extremely unbalanced data) 可 能 会 


F (genome-wide association study, GWAS) 小 


及 百 万 、 干 万 个 位 点 ， 如 果 无 法 严格 控制 一 类 错误 ， 那 么 假 阳 性 位 点 的 绝对 数量 会 很 大 ， 


将 对 后 续 研 究 造成 严重 干扰 和 资源 浪费 。GWAS 极端 不 平衡 数据 主要 体现 在 : (0) 自 变 
低频 甚至 罕见 变异 ， 例 如 次 等 位 基因 
本 库 (UK Biobank, UKB) 级 别 数据 中 次 等 位 基因 
居 ， 但 事件 的 比例 较 低 。 例 如 ， 在 病例 -对 照 研究 中 ， 病 例 与 对 照 比 
， 疾 病 的 发 病 率 较 低 ， 低 于 5%; (3) 因 


Q) 
例 极度 不 平衡 ， 


因 变 量 是 二 分 类 数 


超过 1:20; 


量 是 
频率 (minor allele frequency, MAF) < 0.01 或 英国 生物 标 
计数 (minor allele count, MAC) « 400 的 位 点 ; 


在 队列 下 


FHI 


变量 是 生 


存 数 据 ， 但 删 失 比例 严重 ， 超 过 9099 B 9。 因 此 ， 本 文中 极端 不 平衡 数据 定义 为 自 变量 和 / 
或 因 变量 指标 的 取 值 呈现 严重 比例 失衡 的 数据 。logistic 回归 模型 、Cox 比例 风险 模型 等 参 
数 假设 检验 的 经 典 统计 量 的 应 用 前 提 是 ， 大 样本 条 件 下 ， 统 计量 收敛 于 正 态 分 布 或 卡 方 分 
布 。 然 而 ， 对 于 极端 不 平衡 数据 ， 即 便 在 大 样本 条 件 下 ， 经 典 统 计量 仍然 无 法 有 效 收敛 于 
对 应 分 布 ， 难 以 控制 一 类 错误 5 
本 文 内 容 结构 如 下 : 第 一 节 回 顾 和 常见 统计 模型 参数 假设 检验 的 理论 基础 ; 
经 典 统计 方法 无 法 用 于 极端 不 平衡 数据 统计 分 析 的 理由 ;第 三 节 着 重 介 绍 Firth 校正 和 鞍点 
近似 (saddlepoint approximation, SPA) 方 法 ， 阐 述 两 种 方法 处 理 极 端 不 平衡 数据 的 原理 和 效果 。 
第 四 节 简 介 基 因 组 学 极端 不 平衡 数据 常用 软件 。 
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一 、 经 典 模型 参数 假设 检验 的 理论 方法 回顾 
1 经 典 的 大 样本 检验 
首先 ， 我 们 回顾 三 种 渐 近 等 效 的 大 样本 假设 检验 方法 ， 即 Wald 检验 、 似 然 比 检验 和 得 
分 检验 。 这 三 种 方法 以 极 大 似 然 估 计 和 中 心 极限 定理 为 依据 ， 构 造 大 样本 下 服从 正 态 分 布 
或 卡 方 分 布 的 检验 统计 量 ， 在 参数 模型 的 假设 检验 中 起 核心 作用 。 
(1) Wald 检验 
当 样 本 量 较 大 时 ， 总 体 参数 y 的 极 大 似 然 估计 9 在 适当 的 正则 条 件 下 会 依 分 布 收敛 于 正 


态 分 布 [6 , Bl: 


YEN 
se{ĝĵ} 


式 (1) 中 ，se{ 分 是 极 大 似 然 估计 ?的 标准 误 。 对 于 极 大 似 然 佑 计 ， 相 应 的 se{ 了 } 共 有 性 质 : 


-5N(0,1) (1) 


随 着 样本 量 增加 ， 收 敛 于 1/VIn(Y)。 此 处 的 (yy) 被 称 为 期 望 信息 量 或 Fisher fA RU). Wt 
训 .Xs 来 自 于 总 体 分 布 F(x|y) 的 随机 样本 ， 则 太 () 的 具体 公式 如 下 : 
L() = E( flog (2) 
NM d. 
-- E(grskeb(e)) — (对 于 指数 族 恒 成 立 ) 


=- E(X Patoss(e.9)] 
- -nE (Z zlogf(a1)) 


式 (2) 中 的 L(y|x) = Tf ily) AWA RB log 为 以 e 为 底 的 对 数 。 由 于 真实 的 总 体 参 


BY ETE ARAB, AERO F By Fe RR ETER, AES lsefP} = 
1N). 

HIE mR RH: y = yo， 则 可 构造 式 (3) 中 的 Wald 检验 统计 量 ， 当 吉成 立时 ， 
该 统计 量 在 大 样本 下 服从 标准 正 态 分 布 : 


— Yo 
Zwaa = at 4 } (3) 


以 显著 性 水 平 q 进 行 双 侧 检验 时 ， 若 |Zwara| > Z-a NEA BZA PML. Zaz 
表示 标准 正 态 分 布 的 (1-g/2)X100% 分 位 数 。 同 样 ， 通 过 了 + 1.96sefP KBE By) 95% 置 
信 区 间 。 将 Wald 统计 量 取 平 方 可 得 到 卡 方 统计 量 的 形式 ， 即 aa = ( — 9)? /var(f] = 
(一 Yo)?1n(7?)， 其 服从 自由 度 为 1 的 卡 方 分 布 。 若 X%aq > X-a) EHH. 

(2) 似 然 比 检验 

似 然 比 检验 (likelihood ratio test, LRT) 是 另 一 个 广泛 使 用 的 假设 检验 方法 ， 当 对 仅 有 的 
一 个 总 体 参数 y 进 行 假设 检验 Ho:y = yo 时 ， 似 然 比 统计 量 的 定义 为 : 


由 此 可 见 ， 似 然 比 统计 量 为 当权 成 立时 的 似 然 函数 除 以 最 大 的 似 然 函数 ， 显 然 分 母 中 
y 的 取 值 便 是 极 大 似 然 估计 ?了 。 由 于 似 然 函 数 为 非 负 数 ， 似 然 比 统计 量 的 值 域 为 0 到 1. 25 
yo 时 ， 似 然 比 统计 量 达 到 最 大 值 1， 说 明 有 充足 的 证 据 支 持 Ho; 反之 当 似 然 比 统计 量 接 
0 时 ， 说 明 在 条件 下 出 现 该 样本 的 可 能 性 极 低 ， 倾 向 于 拒绝 Ho。 
在 大 样本 下 ， 对 数 似 然 比 检验 统计 量 收敛 于 自由 度 为 1 的 卡 方 分 布 ， 即 : 


= 


L(yolx) _ 


Wa 2[logZ (4| a) — logL (yol a) ] 一 X2(1) (5) 


Xtrt = — 2log 


给 定 显著 性 水 平 Kc， 当 xzgr > Xall) H HEH. RE EIR AVE ie HI UAE SUY IT] 
$e BiU (profile likelihood) &fzi EX JiS Wii 959656 Bb LIMA Efe E 7491 : 
2|logL L(^| z) —logL (Yupper | 2) | 


= 2[logL (^|) — log L (i12) | (6) 
= xà. 95 — 3.84 


求解 式 (6) 中 的 yupper 和 Yiower 可 得 95% 置 信 区 间 的 上 、 下 限 。 
(3) 得 分 检验 
得 分 检验 (score tesb 又 被 称 为 拉 格 明 日 乘 子 检验 (Lagrange multiplier test)。 定 义 得 分 统计 


量 (score) 为 对 数 似 然 函 数 的 一 阶 偏 导 ， 即 : 


四 | 


S(y)= Frlogb Gl) - 2^5. 7 log fni») (7) 


"AH: y = yo 成 立时 ，S(yo) 在 大 样本 下 服从 均值 为 0， 方 差 为 (yo) 的 正 态 分 布 。 因 此 ， 
得 分 检验 的 检验 统计 量 为 


S (yo) "E S (yo) 
R TO CIC ay ADD ®) 


与 Wald 检 验 类 似 ， 在 大 样本 下 Zscore 服 从 标准 正 态 分 布 。 当 以 显著 性 水 平 g 进 行 双 侧 检 
验 时 ， 若 |Zscore| > Zi1_ay2， 则 拒绝 Ho。 对 Zscore 取 平方 可 得 到 卡 方 统计 量 的 形式 ， 即 
Xscore = S(Yo)?/In(Yo)， 当 Xscore > Xf-a(1) 时 ， 拒 绝 Ho。 

(4) 多 元 统计 中 的 三 种 大 样本 检验 
推广 到 随机 变量 中 含有 多 个 未 知 总 体 参 数 的 场景 。 记 6 = (0,...,0,) 是 p 维 未 知 的 总 
体 参 数 向 量 ， 参 数 的 极 大 似 然 估计 6 在 大 样本 下 会 依 分 布 收敛 于 均值 为 6， 方差 - 协 方差 矩 
阵 为 ,~1(0) 的 多 元 正 态 分 布 中 。L,(9) 此 时 被 称 为 信息 和 矩阵， 具有 如 下 形式 : 


1,(0) = LS veh (81a) E lost (6\x)| B 


" O?logL (0lz) 
= al 06007 ) 0) 


2 ng (orae) 


0000" 


如 果 对 p 维 的 参数 向 量 9 做 双 侧 假设 检验 : Ho: 0 = 0o, H1: 0 4 96o， 三 种 大 样本 检验 统 
计量 的 构造 如 下 : 


(1) Wald 检验 统计 量 : xWala = (8 — Oo)" 1, (8)(0 — 09): 


(2) 似 然 比 检验 统计 量 : Xr 2[log L@|x) — log L(69|x)]: 


(3) 得 分 检验 统计 量 : X&eore = S0) In (9o)S(6o)。 
在 原 假设 下 ， 以 上 三 个 检验 统计 量 均 渐 近 服从 自由 度 为 p 的 卡 方 分 布 。 在 给 定 显著 性 
水 平 wc 下 ， 当 检验 统计 量 大 于 好 Qq)hTBA Je He. 


假设 只 对 p 个 未 知 总 体 参 数 中 某 一 个 分 量 进行 假设 检验 ， 记 8 = (ESR PAR 
Y 


检验 第 p 个 参数 ， 其 余 参 数 B 为 见 余 参数 (nuisance parameter)。 对 y 进 行 如 下 假设 检验 : 
Hg: y = Yo Hy Yo， 则 三 种 大 样本 检验 统计 量 的 构造 如 下 : 


(1) Wald 检验 统计 量 :x&aia = (一 yo)?/In (0), 


(2) 似 然 比 检验 统计 量 : yeep = 2[log LB, flx) — log L(Bo, Yol]: 


(3) 得 分 检验 统计 量 : Hicone = SOo)zm * (8),, 


ERP, 8 = (4) 为 所 有 参数 的 极 大 似 然 估计 值 ， 而 06 = (名 ) 为 Ho:y = REN, I 


余 参 数 B 的 极 大 似 然 估计 值 。1-1(B),, 和 1/14+(Bo),, 为 矩阵 对 角 线 上 第 p 个 元 素 ， 分 别 代 


lm 


表 极 大 似 然 估 计 了 和 得 分 统计 量 5S(yo) 的 方差 分 量 。 以 上 检验 统计 量 大 样本 下 均 服 从 卡 方 分 
布 ， 自 由 度 为 10。 
表 1 汇总 了 三 种 检验 统计 量 的 性 质 。 

表 1. 三 种 渐 近 检验 统计 量 的 性 质 


全 局 检验 Ho:9 = Og 局 部 检验 Ho:y = yo 与 似 然 函数 的 关系 
Y = Yo? n 
Wald 检验 (8 — 00) 1,(8)(8 - 65) 2 x2 £ = 9; 3d 只 需要 优化 无 约束 模型 
" pp 
似 然 比 检验 2[logL(8|x) — logL.C091x)] 3x 2[logL(8|x) — log (8,|x)] 5 需要 同时 优化 有 约束 和 无 约束 模型 
得 分 检验 S(09) 1, 1(00)S(00) B xi S(yoY 1s * (85), 343 只 需要 优化 有 约束 模型 


Ho:0 = 90， 例 如 检验 go = (06,..,0,) = 0; 局 部 检验 是 指 在 控制 了 其 他 协 变量 后 ， 检 验 
某 一 个 参数 y。 三 种 检验 都 和 似 然 函 数 密切 相关 ， 对 所 有 参数 进行 极 大 似 然 估计 的 模型 被 称 
为 无 约束 模型 ， 而 Ho 成 立时 ， 仅 对 宛 余 参 数 进行 极 大 似 然 估 计 的 模型 被 称 为 有 约束 模型 。 
由 此 可 见 ，Wald 只 需 拟 合 无 约束 模型 ， 得 分 检验 只 需 拟 合 有 约束 模型 ， 而 似 然 比 检验 需 同 
时 拟 合 两 个 模型 。 


1.2 logistic 回归 模型 


归 模型 如 下 ; 


GWAS 研究 中 常用 的 logistic 


logit [Pr (Y;=1|X;,,G;)]= logit m; = X; B+GiY, t=1,2,.…,n (10) 


X10), Y, = 1 或 0， 分 别 代 表 第 i 个 个 体 是 病例 或 对 照 ，zti 表 示 该 样本 是 病例 的 概率 ， 
Xi 是 k x 1 维 截 距 项 与 协 变 量 ，Gi 代 表 第 i 个 个 体 待 检验 的 遗传 位 点 ， 编 码 为 次 等 位 基因 数 
H. G; = 0,1,2。B 是 k x1 维 协 变 量 的 系数 向 量 ,，y 是 位 点 的 效应 系数 ， 等 于 其 优势 比 
(odds ratio, OR) 的 对 数 。GWAS 关注 该 遗传 位 点 是 否 与 疾病 相关 ， 则 原 假 设 Ho:y = 0。 


H 


归 系数 p 和 y 均 可 采用 极 大 似 然 方法 进行 估计 。 记 9 = (人 为 kt1 维 的 系数 向 量 


y,X,G 分 别 代表 了 每 一 行为 yj, XG nT ÉE E, logistic 回归 的 似 然 函数 如 下 : 


代 。 


Ti 


In 


L(61y,X.G)- | [Pr (x; 2 16,x.6)- | [x*(1-)* (11) 
i=1 i=1 


XT BULA PR — MEA eo SH MILI, ARR EI TET: 


“ à 
s(0)— 2.3108 Pr (Y; = y;|X;,8,G:,7) 
i=1 


E 
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(12) 


当 画 成 立时 ， 零 模型 的 到 可 用 极 大 似 然 估计 人 多。 = exp(X?Bo)/[1 + exp(X1 Bo) ttr & 
由 此 可 见 ，R 为 零 模 型 的 残 差 向 量 R = (Y, — figi... Yn — Rony o 
寺 数 似 然 函 数 的 二 阶 偏 导 为 Fisher 信息 矩阵 ， 即 
(3 (01y,X, 2) 
0000" 


A 0? log Pr(Y; = y;|0,X;,G,) 
sp 806807 
Sues did 

 1G"WX G™WG 


>4 


I,(0)— 


(13) 


ROPP IW ALY, KIR- BU Dm (A m APER TCR IIO PRG, IR 


TIEN BECK ADLER f Te; FS TK S XR DR ES SET DAT BO HY Fisher 信息 矩阵 的 逆 矩 阵 


1(9)， 其 对 角 线 上 的 第 Kt1 个 元 素 ( 即 参数 y 对 应 元 素 ) In (Oke = 1/[GTWG — 


G'WX(X'Wx) !x'WaG]. 


结合 以 上 推导 ，logistic 回归 的 三 种 大 样本 检验 统计 分 别 由 下 式 给 出 : 
(D Wald 检验 统计 量 ; 


Zwaia = x m i 
T PRIORE (14) 


-4ya"Wa-G"Wx(x"WXx) x'Wa 


Q) 似 然 比 检验 统计 量 ; 


Xx = 2|logL (8,4 


y, X,G) — log! (Bo, yo = 0 |y, X. G)] 


= aS DERE 十 log(1 一 | — Yes ü E + log(1 iw} (15) 


ot (1 — io) , =a 
=2) [utes a GE | log Ti, 


(3) 得 分 检验 统计 量 : 


$(0) —— $(0) 

Vear(S(} fin, (à; 
_ GTR 

VG"W.G—-G"W,X(X"W,x) X'W.G 


k+1,k+1 


(16) 


值得 注意 的 是 ， 在 GWAS 无 缺失 数据 的 情况 下 ， 得 分 检验 只 需 拟 合 一 次 固定 协 变量 组 
合 的 空 模 型 ， 便 可 得 到 式 (16) 中 的 重要 参数 R 和 Wo， 故 运算 效率 很 高 ， 在 大 规模 GWAS 中 


£I iz JU Um, dd = G —X(XTWX) XTWoG， 为 协 变量 调整 的 位 点 向 量 ， 则 可 


计算 另 一 种 形式 的 得 分 统计 量 ( 式 17) 及 其 渐 近 方差 ( 式 18)。 两 种 方法 所 得 检验 统计 量 Zscore 
完全 相等 1， 但 是 具有 更 为 简洁 的 形式 。 


5(0)= YDG: Y,- ão) -G7R a7 
var(S(0)) = GTW,G (18) 
Fina OM) an CIE (19) 


Vvar{S(0)}  VGTW,G 
因此 ，GWAS 中 logistic 回归 一 般 采 用 如 下 步骤 进行 得 分 检验 ， 以 提高 计算 和 存储 的 效 
率 。 首 先 ， 拟 合 结局 与 协 变量 的 logistic 回归 模型 ， 算 得 如 R、Wo; 然后 ， 拟 合 待 检验 位 点 
和 协 变 量 的 线性 回归 模型 ， 算 得 G6; 最 后 代入 式 (19)， 便 可 得 Zscore。 若 分 析 不 同位 点 ， 只 
需 更 新 G 即 可 ， 从 而 避免 反复 进行 极 大 似 然 法 估计 参数 ， 计 算 速 度 较 快 。 但 是 ， 该 方法 不 
适用 于 构造 置信 区 间 。 


HE 


1.3 Cox 比例 风险 模型 
Cox 比例 风险 模型 可 用 式 (20) 表 达 : 


此 处 ，4(t; Xj, Gi) 表 示 样 本 在 t 时 刻 的 风险 函数 ，X0(t) 表 示 基 线 风 险 函 数 。 风 险 函数 
4() 代 表 样 本 在 t 时 刻 发 生 结 局 事件 的 瞬时 概率 。 同 样 ， 原 假设 为 Ho:y = 0。 该 模型 的 似 然 
函数 由 于 包含 未 知 的 风险 函数 416， 无 法 直接 进行 极 大 似 然 估 计 。 因 此 ，Cox 提出 采用 偏 似 
然 函数 估计 目标 参数 3。 令 了 R(t) = :7 三 如 表示 在 t 时 刻 承受 风险 的 样本 集合 ， 则 似 然 函 
数 中 40 可 以 消 掉 ， 得 到 偏 似 然 函 数 : 


A 


£ exp(X7 B 4- G; 
L, (8.1 A,T.X,G)- [| PLO 
1] 2, exp(X? 8-- Gy) 


jeR(t) 


(21) 


观测 到 的 生存 结局 由 两 部 分 组 成 : COE SE PUA, ARIA; = 1 或 0; @ 
生存 事件 T;。 式 (21) 仅 取决 于 事件 时 间 的 秩 次 ， 且 没有 直接 使 用 删 失 和 未 删 失 的 事件 时 间 数 
值 ， 因 此 被 称 为 偏 似 然 函 数 (partial likelihood fonction)。 利 用 偏 似 然 函数 可 对 模型 中 参数 有 
和 7 进行 极 大 偏 似 然 估 计 。 当 生存 数据 中 存在 较 多 打 结 (ties) 数 据 时 ， 即 存在 大 量 相 同 取 值 
的 生存 或 删 失 数据 ， 则 可 采用 Breslow's 近似 或 Efron's 近似 等 方法 校正 偏 似 然 函 数 04。 

偏 似 然 函 数 的 y 一 阶 偏 导 为 得 分 统计 量 。 当 Ho 成 立时 ， 得 分 统计 量 为 : 
OlogL, (B, y) 

a 


S(0)= 


5 G;exp(X;' B) 


= = A; G, JER(t) 
2. Y exo (Xj B) 


JER) 


= YDG [A - H(,X.8)] 


~ DDG [A - n(x.À)] 


E SGR, -G"R 
i=l 


(22) 


此 处 ，H(ti,X,pB) 是 累积 风险 函数 ，Ri = A; H(t X, Bo) ASR M PHAM MRE 


(martingale residual). FJ WL, Cox 回归 模型 与 logistic 回归 模型 的 得 分 统计 量具 有 一 致 的 形式 


[5], 


对 数 偏 似 然 函数 二 阶 偏 导 为 Fisher 信息 矩阵 。 记 9 = (5 ) 为 tt1 维系 数 向 量 ， 并 计算 仿 


JU EA PR BAY) Fisher fei BER: 


2 
1,(6)=5| Ó a ord 


00007 
2 pad See] 
 (GTQX G"QG 


此 处 ，Q@ 是 结局 事件 的 方差 - 协 方差 矩阵 ， 有 具体 公式 为 : 


Q3) 


Q= Y; ys [o (dag (Ad — AA] 
A; — 150 * exp(XB + Gy) 
w; (0) = 124, exp (XB + Gy) = 3 5 A; 


(24) 


式 中 ， Ine 为 一 个 n 维 的 示 性 向 量 ， 用 于 指示 生存 时 间 长 于 第 i 个 个 体 的 样本 j， 即 


U:T; 2 tid7g 1, RZWA 0. diag{- AVA iA) EAD ARTO RIED fü FEE. AIA SY 
R, EXAMEN ICR BUSTA, thw; (8) ASA, FE PICK ZAM 

于 是 ， 三 大 渐 近 检验 统计 量 的 形式 如 下 : 

(D Wald 检验 统计 量 : 


Zwaid = m = d 
A VT (8), (25) 


-4Ja*Qa -a*Qx(x"Qx) Xroc 


(2) 似 然 比 检验 统计 量 : 


A,T,X,G) — logL, (boyo =0|A,T,X,G) | 


Xer = 2|logL, (8,5 


n | Xo exp( X7 B +G:4) 
= A; XI B4 G^ T By — log| 255? - 
Cale S srt 26 
jER(t) 


w,(6) 


=2) A; 号 B 4- G,5 — XT B, el 


(3) 得 分 检验 统计 量 ( 令 G = G-X(X'QyX) X70@06): 


S(0)  . S(0) 
VivertS(0)E AEB) os 
GTR 


score 


| Ve'àe-erQ.x(x'Q.x) x'à Tm 


GTR 
GTG 

13 混合 效应 模型 

种 群 结构 和 亲缘 关系 是 影响 GWAS 结果 的 重要 混杂 因素 。 混 合 效应 模型 是 常用 的 解决 
FRU S 1618, EF Cox 混合 效应 模型 也 拥有 相似 的 形式 ， 本 文 以 logistic 混合 效应 模型 作 
为 示例 ， 具 体 模型 如 下 


logit [Pr (Y; =1|X,,G,)]=logit v; = Xi 8 - G;y +é, i=1,2,--,n (28) 


与 传统 logistic 回归 模型 不 同 ， 式 (28) 增 加 了 用 于 刻画 种 群 结构 和 亲缘 关系 的 &。€ = 
(5, tn》 为 随机 效应 ， 服 从 多 元 正 态 分 布 Ni(0,TV)。 其 中 , VAn xn 的 遗传 关系 矩阵 


(genomic relationship matrix, GRM)，Tt 是 用 于 刻画 方差 分 量 的 尺度 参数 。 该 模型 中 的 Bp、y 和 


z 三 个 参数 可 使 用 惩罚 准 似 然 (penalized qusai-maximum likelihood, PQML)09] 或 平均 信息 限制 


性 最 大 似 然 (average information restricted maximum likelihood, AI-REML) 算 法 


= Wiki GiRi = Xt 


在 如 条 件 下 ， 得 分 


var{S(0)} = GEG — G'$5X(XT£oX) X'£gG. fJ Ho AK EF ft BH 


L^ T wl = 


统计 量 


为 : S(0) 


[diag{foi/(1 — 9) ! tV] * 


异 和 GRM。 同 样 ， 得 分 检验 统计 量 Zscore 收 敛 到 标准 正 态 


对 于 生存 分 析 的 混合 效应 模 


Wo， 便 可 进行 得 分 检验 。 


二 、 极 端 不 平衡 数据 对 经 典 统计 量 分 布 的 影响 


通过 模拟 试验 结果 阐 


I 述 极端 不 平衡 数 


。 由 此 可 见 ， 数 


分 布 。 


Gi 


(Y; — îoi) > 


Von AT iit 


其 方差 为 : 


FEX. 3 = 


Fi IUS SERO ELT ARE 


型 ， 将 R; TI (Y; m ftoi) EHE, H 


日 前 文中 的 Q@o 蔡 代 


居 对 经 典 统 计量 


不 考虑 协 变量 ， 假 定位 点 和 结局 存在 如 下 关系 : 


分布 的 影响 。 简 单 


logit [Pr (Y; =1|X,,G;)]=a+G;x B, i=1,2,---,n 


起 见 ， 模 型 暂 


Q9) 


参考 UKB 中 的 肺癌 数据 设置 模拟 情景 ， 即 假设 300,000 人 的 队列 中 肺癌 基线 患 病 率 为 


1%, 23 3,000 例 病 例 ， 故 截 


分 别 设置 两 种 情景 : 


位 点 G6 的 MAF 分 别 为 0.01、 


0 生成 模拟 样本 ， 分 别 月 


距 项 a = 一 4.595。 


待 检验 位 点 G 采 | 


模拟 试验 重复 100 万 次 ， 获 得 各 个 统计 量 的 经 验 分 布 ， 并 阴影 


如 图 1 所 示 ， 当 MAF 取 0.01 时 ， 检 验 统 计量 
较为 接近 。 双 侧 检验 尸 值 的 QQ 图 


1。 如 图 2 所 示 ， 当 MAF 下 降 到 0.001 时 ， 


分 布 严 重 偏离 标准 正 


数 极 少 ， 检 验 统 计量 


三 也 呈现 


可 能 存在 严重 的 假 阳性 。 


态 分 布 或 y? f 
出 明显 的 离散 化 趋势 。 此 外 ， 基 因 组 膨胀 因子 大 于 1， 也 表明 


由 于 自 变 量 和 相 


EX 


0.001， 分 别 代 表 低 频 、 
H Wald 检验 、 似 然 比 检验 和 得 分 检验 对 位 点 G 的 系数 进行 假设 检验 。 


EN 


ZU A Hemp E. 


Za 


E 
里 


量 的 经 验 分 布 与 对 应 的 标准 正 态 
显示 ， 基 因 组 膨胀 因子 (genomic inflation factor) 4 趋 近 于 


即便 样本 量 已 达到 30 万 人 ， 检 验 统计 量 的 经 验 


相 加 模型 编码 ，G; = 0,1,2. 
罕见 变异 。 令 B= 


下 2.5% 分 位 数 。 
分 布 或 y? A 


的 取 值 为 离散 型 变量 且 事 件 


(A1) Wald 检 验 统计 量 分 布 (A2) Wald 检 验 QQ 图 


ec 
EN 一 标准 正 态 分 布 g 
jj 一 经 验 分 布 E 
oe [I 
E e 
o 
[7] 
ô 
-1.96-1.89 0 1.962.01 
Expected —logi9(p) 
(B1) 似 然 比 检验 统计 量 分 布 (B2) 似 然 比 检验 QQ 图 
ec 
eo 
H 一 di 8 
Hi 一 sub E 
= 2 
Q 
ri 
O 
0.2 3.84.85 0 2 4 6 
Expected - logio(p ) 
(C1) 得 分 检验 统计 量 分 布 (C2) 得 分 检验 QQ 图 
T ean. 6 - 
Q 
© 
R 一 标准 正 态 分 布 9 
i 一 经 验 分 布 
a 经 验 分 布 E 
E 2 
O 
a 
ô 
-1.96-1.90 0 1.96 2.02 0 2 4 6 


Expected - logio(p ) 


图 1. 轻 度 不 平衡 数据 三 种 检验 统计 量 的 经 验 分 布 


(情境 1: 样本 量 n = 300,000， 患 病 率 r = 1%， 位 点 M4F = 1%) 


(A1) Wald 检 验 统计 量 分 布 (A2) Wald 检 验 QQ 图 


€ 
ps 一 标准 正 态 分 布 g 
x 一 steht E 
ae [I 
E 2 
A 
Q 
O 
-1.96 -1.61 0 1.96 2.10 
Expected — logi9(p) 
(B1) 似 然 比 检验 统计 量 分 布 (B2) 似 然 比 检验 QQ 图 
€ 
e 
aX — x^ Afi 3 
一 经 验 分 布 i 
oO 
2 
oO 
o 
Q 
oO 
0.2 3.88.84 

Expected - logio(p ) 

(C1) 得 分 检验 统计 量 分 布 (C2) 得 分 检验 QQ 图 
€ 
$ 
Ex 一 标准 正 态 分 布 9 
i 一 经 验 分 布 i 
ea o 
Es 2 
Oo 
[7] 
Q 
oO 

-1.96 -1.70 0 1.96 2.14 0 2 4 6 


Expected - logio(p ) 


图 2. 极度 不 平衡 数据 三 种 检验 统计 量 的 经 验 分 布 
(情境 2: 样本 量 n = 300,000, BAET = 1%， 位 点 M4F = 0.1%) 

如 果 错 误 地 以 标准 正 态 分 布 或 y? 分 布 的 界 值 作为 假设 检验 统计 量 的 界 值 ， 那 么 三 种 检 
验 无 法 严格 控制 一 类 错误 ， 如 表 2 所 示 。 一 类 错误 发 生 膨胀 的 程度 与 检验 水 准 的 取 值 相关 。 
当 检 验 水 准 设 为 0.05 时 ， 在 30 万 人 的 样本 下 ，MAF 无 论 为 0.01 还 是 0.001， 双 侧 检验 的 第 
一 类 错误 尚且 接近 0.05。 然 而 在 GWAS 研究 中 ， 因 为 需要 分 析 多 个 位 点 ， 所 以 进行 一 类 错 
误 多 重 校正 。 以 分 析 1000 个 位 点 为 例 ，Bonferroni 方 法 校正 后 的 检验 水 准 设 为 5 x 10-5 时 ， 
此 时 在 MAF=0.01 的 情景 下 ，Wald 和 得 分 检验 一 类 错误 分 别 膨胀 至 6.60x105 (132%) 和 
7.22x10? (144%)， 似 然 比 检验 的 一 类 错误 控制 较 好 ， 而 在 MAF=0.001 的 情景 下 ，Wald 和 
得 分 检验 一 类 错误 更 是 膨胀 至 2.23 104 (446%) All 7.22 10° (702%)， 似 然 比 检验 一 类 错误 控 
出 过 于 严格 ， 偏 保守 。 若 位 点 MAF 更 小 或 基线 患 病 率 更 低 ， 上 述 现象 更 明显 下 和。 主要 
原因 是 ，GWAS 试图 挑选 P 值 特别 小 的 位 点 ， 极 端 不 平衡 数据 分 析 所 得 三 种 统计 量 的 经 验 

分 布 与 标准 正 态 分 布 或 ?分 布 在 分 布 尾 部 的 差异 相对 而 言 较 大 ， 这 将 造成 严重 的 一 类 错误 


a 


膨胀 或 一 类 错误 保守 导致 的 检验 效能 低下 。 此 外 ， 由 于 极度 不 平衡 数据 的 检验 统计 量具 有 
偏 态 的 特点 ， 因 此 单 侧 检 验 的 一 类 错误 呈现 更 为 严重 的 偏 倚 。 
He 2. 极端 不 平衡 数据 的 三 种 方法 假设 检验 第 一 类 错误 率 


双 侧 检验 左 侧 检验 右 侧 检验 
显著 性 水 平 a MAF 方法 (显著 性 水 平 a) (显著 性 水 平 w/2) (显著 性 水 平 a/2) 

拒绝 域 第 一 类 错误 拒绝 域 第 一 类 错误 拒绝 域 第 一 类 错误 
0.01 Wald 检验 Zwaal > 1.96 488x107 Zwad < —1.96 2.08 x 1074 Zwaa > 196 281x107 

0.01 AA EC US Aber > 3.84 5.04 x 107? — — — — 
0.01 得 分 检验 Zscore| > 1.96 4.95 x 107? Zscore < —1.96 2.12 x 107? Zscore > 196 — 293x107 
0.001 Wald 检验 Zwaal > 1.96 3.34x107 Zwaa < —1.96 «1.00 x 1079 Zwaa? 196 — 334x107 

0.001 AR LR Arr > 3.84 4.95 x 107? — — — — 
0.001 得 分 检验 Zscore| > 1.96 5.11 x 107? Zscore < 1.96 1.59 x 107? Zscore > 1.96 — 352x107? 
5x 1075 0.01 Wald 检验 Zwaal>4.06 — 6.60x1075  Zwaa < —4.06 4.00 x 1075 Zwaa > 406 6.20 x 10-4 

0.01 AA EC US. Xirr > 16.45 5.30 x 1075 = = = = 
5x 1075 0.01 得 分 检验 Zscorel > 406 7.20x10-5 Zscore < —4.06 4.00 x 1075 Zscore > 406 6.80 x 1074 
0.01 Wald ff en Zwaal > 406 2.231074  Zwaa < —4.06 «1.00 x 1079 Zwaia > 406 — 223x107 

0.001 AA EC Xirr > 16.45 2.70 x 1075 一 a 一 — 
5 x 1075 0.001 得 得 分 检验 Zecore| > 4.06 3.51x104 Zoue<-406 — «1.00x 1075 Zecore > 4.06 — 351x107 


=. GWAS 研究 中 极端 不 平衡 数据 的 常用 校正 方法 

3.1 Firth 校正 方法 

虽然 在 大 样本 下 ， 参 数 的 极 大 似 然 估计 6 会 收敛 于 真 值 9， 但 是 当 样本 量 较 小 或 数据 极 
端 不 平衡 时 ，6 作 为 9 的 估计 往往 是 有 偏 的 ， 偏 倚 程 度 取决 于 得 分 函数 5(9) 的 曲率 中。 为 了 
消除 0 Qn D) Er. Firth 提出 采用 惩罚 似 然 函 数 ( 式 30) 进 行 参 数 估计 ， 同 时 试图 控制 一 类 


错误 ?2 


‘SI 


N 


L* (0) - L(0) x|L, (0)|°” (30) 


RP, (UR Be C A TTE, (0) 05 Je rs OBAT UIE, BER Jeffreys 先 验 。 
对 惩罚 似 然 函 数 求 一 阶 偏 导 可 得 Firth 校正 的 得 分 统计 量 为 : 


S*(0) = 8() + 5tr(1, (6) [ƏL (0)/00,]) (31) 


em 


Qx 指 的 是 参数 向 量 9 中 的 第 k 个 分 量 ， 也 即 前 文中 的 y。 式 (10) 中 的 logistic 回归 模型 ， 
相应 的 Firth 校正 得 分 统计 量 为 : 


SOS OE Q2 — 0) (32) 


web, njélBTABEEH = WWM2X(XTWX)-1XTW1/2 的 第 i 个 对 角 元 素 ，W 是 8 的 方差 - 协 
方差 矩阵 。 从 公式 (32) 可 见 ，Firth 校正 同时 从 因 变 量 和 自 变量 两 个 角度 出 发 校正 得 分 统计 
EE。 从 因 变 量 的 角度 ， 当 zo 越 偏离 0.5, Firth 校正 力度 越 大 ;从 自 变 量 的 角度 ，hi 反 映 自 变 
量 的 第 i 个 观测 值 与 自 变量 平均 观测 值 间 的 标 化 距离 ， 最 远 为 1， 最 近 为 0。 当 hi 越 接近 1 


lim] 


7 


Firth 校正 力度 越 大 。 当 因 变 量 和 自 变 量 同时 为 极端 不 平衡 数据 时 ，Firth 法 同时 从 两 


校正 ， 校 正 力度 较 大 。 
(1) Firth 校正 的 Wald 检验 和 得 分 检验 


极 大 化 式 (30) 可 得 Firth 校正 的 参数 估计 6*。 人 的 方差 仍 用 11( 信 ) 的 对 角 线 元 素 进 行 


估计 ， 因 为 其 是 [-0?logL*/(9B)*] 的 一 阶 近似 。 据 此 ， 可 算得 Wald 检验 统计 量 ， 


个 角度 


同 理 也 


可 以 用 Firth 校 正 得 分 统计 量 和 I,,“1(6。 ) 进 行 得 分 检验 。 但 是 由 于 数据 极端 不 平衡 的 情景 


统计 时 呈现 非 对 称 分 布 ，Firth BEL 


[25], 


(2) Firth 校正 的 似 然 比 检验 


Firth 校正 最 常 


in| 


XPLRT == 2|logL* (8*,4* YX, 
=2|logL (8*,5* |y, X, G) + loglI, (8,3155... — log L (Bt, vo =0|y, X, G) — log, (B)|° 5. | 


上 式 


， 疏 和信 是 无 约束 模型 的 极 大 惩罚 似 然 


以 然 函 数 蔡 换 为 惩罚 似 然 ， 即 : 


G) —logL" (Bi, o — 0| y, X,G)] 


ER] Wald 检验 以 及 得 分 检验 仍 无 法 有 效 控 制 第 一 类 错误 


j 的 惩罚 似 然 比 检验 ， 该 方法 也 可 构造 轮廓 似 然 置信 区 间 。 该 检验 统计 
量 是 将 似 然 比 检验 中 的 1 


G3) 


"i. MOBS AAS A A AE TI 


然 估 计 。 相 似 地 ， 轮 廓 似 然 的 95% 28 fea XX Tal AY ERE ct OR AlogL* (Bo, Yo) = logl*(B*, $^) 一 


3.84/25 


的 yo 得 到 。 P 


F Firth 校正 的 似 然 比 检验 需要 通过 和 迭代 得 到 有 约束 模型 和 无 约束 模 


型 的 惩罚 似 然 函 数 极 大 值 ， 故 运算 速度 较 慢 。 


3.2 


鞍点 近似 方法 


鞍点 近似 方法 通过 数据 从 提取 分 布 的 信息 ， 构 造 检验 统计 量 的 经 验 分 布 ， 从 而 根据 Ho 


假定 下 的 经 验 分 布 构造 P 值 。 定义 M(t) = E(e™) = jee FG) dx 7355 RE ER Zt (moment 


generating function, MGF), K(t) = log(M(t))7J £&TA5& BEBÉ Zt (cumulative generating function, 


CGF). Ate At REA oo p 


总 体 均 值 和 方差 : 


(AP; 


XT ES EE AEE PB BL AE f, EL ZR AR n BR ICT ) 


E 质 ， 即 累积 量 母 函数 在 0 位 置 的 一 阶 、 二 阶 导数 分 
EE 
E(X)— K'(0)— 8t Lue 
var(X) — K"(0) — 9 


} 别 对 应 


(34) 


鞍点 近似 方法 实现 高 精度 近 


rae | ü ^ 
olo = tog(4)} for c#4 


F(2)=P(X <2) (35) 


1 K" (0) 
2 6/20 K"(0)” 


其 中 , W= sgn(é) l2(ix -KG aifi ©, son 2S MH, MEEK Ê) = 
x 的 解 。 正 态 近 似 相当 于 只 用 了 前 两 个 累积 量 : 均值 与 方差 。 鞍点 近似 则 使 用 了 整个 累积 
量 母 函数 的 信息 ， 其 相对 误差 界限 为 0(n-32)， 具 有 更 高 的 精度 。 

对 于 检验 统计 量 9， 双 侧 检验 的 P 值 是 Ho 假定 下 比 当前 统计 量 s 更 极端 统计 量 的 比例 ， 
即 P = P(S > |sD + P(S < -IsD = 1— F(IsD 十 F(-|s))。 极 端 不 平衡 数据 的 统计 量 往往 
非 正 态 分 布 。 因 此 ， 基 于 鞍点 近似 方法 所 得 累积 分 布 函数 FC) 所 得 P 值 到 比 正 态 近 似 法 更 
加 精确 。 

三 大 检验 统计 量 中 得 分 检验 统计 量 的 形式 最 为 简单 ， 因 此 也 最 适用 于 提供 鞍点 近似 构 
造 近似 分 布 。 以 logistic 回归 模型 的 色 成 立时 的 得 分 统计 量 $(0) = DL, Gi(; — mi) = 
?1 Gi(Yi 一 倪 01) 为 例 ， 式 中 只 有 是 随机 变量 ， 服 从 伯 努 利 分 布 (n = 1 的 二 项 分 布 ) 的 ， 即 
得 分 统计 量 9(0) 是 蕊 的 线性 组 合 。 大 随机 变量 X 服 从 二 项 分 布 Bi(Cu mw), HRA T IRIE RER 
数 Mx(t) = (re 二 1 一 mi ， 而 独立 随机 变量 线性 组 合 的 和 矩 母 函数 具有 以 下 公式 ; 


for z—ü 


m 


í (36) 


My = | a = | [e Mx.(a.t) 
i=1 


结合 以 上 公式 得 到 得 分 统计 量 5(0) 的 累积 量 母 函数 及 其 一 、 二 阶 导数 。 此 处 将 x 蔡 换 为 
极 大 似 然 估计 看 即 可 得 到 所 有 logistic 回归 考点 近似 方法 需要 的 统计 量 。 
当 结 局 事件 的 具体 分 布 未 知 时 ， 可 以 借助 经 验 分 布 求解 鞍点 近似 中 所 需 的 统计 量 。 记 


My (t) = Dfl etti Jy AS Ap RES ML, 通过 求 导 运 算 可 估计 经 验 累积 量 母 函数 PI， 见 公式 


(37). 


Kx(t) =logM,(t) = (4 oe 


gu) = Ma) _ 2 d " 
Mx(t) Stet 
goo) = Mes Malt) = [Meee _ Dave 2 e^ "Do me^] 
[Mx(t)P [5 ee 


图 3 比较 了 得 分 统计 量 的 正 态 近似 迁 点 近似 累积 分 布 函数 ， 可 见 鞍 点 近似 方法 可 以 
更 好 地 拟 合 分 布 的 头 尾部 分 。 但 是 ， 正 态 近 似 方法 在 x = 0《〈 相 当 于 鞍点 近似 方法 x = 5 
附近 的 拟 合 性 能 并 不 差 。 鉴 于 正 态 近似 方法 运算 效率 远 高 于 鞍点 近似 方法 ， 因 此 一 般 推荐 
先 采 用 正 态 近似 法 计算 P 值 ， 再 对 特别 小 的 P 值 采用 鞍点 近似 法 获得 更 精确 的 数值 。 除 计 
算 速 度 慢 外 ， 鞍 点 近似 方法 不 适合 估计 效应 系数 7 及 其 标准 误 。 


经 验 得 分 检验 统计 量 分 布 
得 分 检验 统计 量 分 布 ( 正 态 近似 ) 
得 分 检验 统计 量 分 布 (鞍点 近似 ) 


0 
得 分 检验 统计 量 
1.00 


0.99 


-3.0 -2.5 -2.0 -1.5 -1.C 1 2 3 
得 分 检验 统计 量 得 分 检验 统计 量 


图 3. 得 分 统计 量 正 态 近 似 和 鞍点 近似 累积 分 布 函数 的 比较 


3.3 总 结 

本 节 简 述 了 Firth 校正 方法 和 遂 点 近似 方法 控制 一 类 错误 的 核心 思想 : 前 者 通过 校正 极 
大 似 然 估计 的 偏 傈 从 而 提高 参数 估计 的 精度 和 假设 检验 的 性 能 ， 后 者 则 通过 数据 构造 检验 
统计 量 的 经 验 分 布 以 进行 假设 检验 。 两 种 法 均 以 计算 复杂 度 为 代价 优化 假设 检验 的 统计 学 


性 质 。 因 此 ， 一 般 推荐 分 析 极 端 不 平衡 数据 或 经 典 检验 P 值 较 小 可 能 影响 一 类 错误 时 使 用 


避免 GWAS 额外 增加 计算 负担 。 以 模拟 试验 情景 2 为 例 ， 仅 对 P 值 <0.01 的 位 点 分 别 实施 


Firth 校正 方法 和 毅 点 近似 方法 。 如 图 4 中 所 示 ，Firth 校正 的 Wald 得 分 检验 和 得 分 检验 仍然 


显 保守 ， 这 可 能 是 Firth 校正 未 能 在 GWAS 


出 现 了 一 类 错误 膨胀 ，Firth 校正 的 似 然 比 检 验 稍 
方法 中 广泛 应 用 的 原因 之 一 。 而 应 用 最 为 广泛 的 鞍点 近似 校正 的 得 分 检验 能 够 较 好 地 控制 
一 类 错误 。 此 外 ， 得 益 于 罕见 变异 基因 型 的 编码 大 部 分 为 0， 这 部 分 样本 信息 对 得 分 统计 


量 无 贡献 ， 通 过 优化 算法 还 可 以 进一步 提高 得 分 检验 的 计算 速度 。 精 确 检 验 、logF 校正 、 
bootstrap 自助 法 等 其 他 校正 检验 存在 计算 效率 低下 的 缺陷 ， 在 GWAS 中 应 用 更 少 ， 此 处 不 


作 介 绍 ， 详 见 相 关 文 献 C&%' 291, 


o 


(A) Wald 检 验 (Firth 校 正 )QQ 图 (B) 似 然 比 检 验 (Firth 校 正 )QQ 图 
6 Pd . 
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A 

Observed - log4o(p) 
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N 
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(C) 得 分 检验 (Firth 校 正 )QQ 图 (D) 得 分 检验 (鞍点 近似 )QQ 图 
g . 
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op 
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AR 

Observed - log4o(p) 
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图 4. Firth 校正 和 鞍点 近似 方法 的 QQ 


四 、 基 因 组 学 极端 不 平衡 数据 的 常用 软件 简介 


ASI GWAS 中 常用 的 极端 不 平衡 数据 分 析 软 件 ( 表 3) 0 


表 3. GWAS 常用 极端 不 平衡 数据 分 析 软 件 


软件 结局 类 型 尸 值 校正 方法 亲缘 关系 校正 
R 软件 包 logistf 二 分 类 结局 Firth x 
R 软件 包 coxhf 生存 结局 Firth x 
PLINK2 二 分 类 结局 Firth x 
R 软件 包 SPAtest 二 分 类 结局 i eI UL X 
R 软件 包 SPACox 生存 结局 鞍点 近似 x 
SAIGE 二 分 类 结局 鞍点 近似 y 
GATE 生存 结局 鞍点 近似 y 
REGENIE 连续 性 、 二 分 类 结局 Firth、 鞍 点 近似 Y 


4.1 R 软件 包 logistf 和 coxphf。R 软件 中 的 logistf 和 coxphf 分 别提 供 Firth 校正 的 logistic 


H 


回归 和 Cox 


归 ， 并 采用 Wald 统计 量 和 惩罚 似 然 比 统计 量 实现 假设 检验 。 虽 然 Firth 方法 


校正 了 回归 模型 


的 系数 估计 ， 但 使 得 预测 概率 向 0.5 ffir. AE Iogistf 包 中 提供 了 FLAC 和 


FLIC 两 种 方法 来 改善 预测 结果 B0。 
4.2 PLINK2。PLINK2 是 知名 GWAS 数据 分 析 与 处 理 软件 的 最 新 版 本 B0。 该 软件 默认 


处 理 策 略 是 : 9 


回归 。 该 软件 不 能 


E 对 所 有 位 点 进行 logistic 回 归 ， 再 对 模型 不 收敛 的 位 点 采用 Firth 校 正 logistic 


自动 鉴别 极端 不 平衡 数据 ， 推 荐 对 MAC < 400 的 位 点 实施 Firth 校正 外。 


4.3 SPAtest。Dey 等 [将 鞍点 近似 方法 引入 二 分 类 结局 的 GWAS， 并 发 布 了 R 软件 
SPAtest 包 。 该 包 对 得 分 统计 量 位 于 均值 土 2 倍 标准 差 范围 以 内 的 检验 统计 量 采 用 正 态 近 似 ， 
范围 以 外 的 检验 统计 量 采 用 鞍点 近似 ， 以 得 到 假设 检验 的 已 值 。 些 外， 用户 可 以 依据 


Berry-Esseen 7E PEAY IN 


值 进行 判别 。 研 究 者 还 提供 了 FastSPA 的 版 本 ， 针 对 基因 型 中 大 量 的 


0 数据 进行 了 入 


4.4 SPACox. 


法 优化 ， 以 提高 运算 效率 。 
该 方法 是 鞍点 近似 校正 的 Cox 比例 风险 模型 ， 适 合 分 析 大 规模 生物 标本 
库 级 别 的 生存 数据 PH。 


| 


对 于 二 分 类 结局 ， 得 分 统计 量 服 从 二 项 分 布 ， 易 于 构造 累积 量 母 函 


数 。 对 于 生存 数据 ， 蒜 残 差 较 难 直接 得 到 分 布 函数 。 因 此 ，SPACox 直接 估计 经 验 累 积 量 母 


函数 ， 从 而 构造 


S 


量 的 经 验 分 布 以 进行 假设 检验 。 此 算法 可 通过 及 软件 SP4Cox 包 实现 。 


= 


4.5 SAIGE. SAIGE 实现 了 鞍点 近似 校正 的 广义 混合 效应 模型 ， 相 较 于 SP4test 包 中 的 


FastSPA 方法 ， 该 方法 利用 GRM 额外 校正 了 种 群 结构 和 亲缘 关系 叫 。 在 运算 方面 ， 得 分 统 


计量 的 方差 vars(0) = G'£,6 — G'E X (XTE X) X730G = GILG KHIM. Aue, 


SAIGE 首先 用 部 分 位 点 估计 变异 的 比例 ?= oT. ELST MIE SOIT, Rn 
通过 正 态 近 似 或 鞍点 近似 计算 己 值 ， 以 同时 优化 计算 速度 和 统计 学 性 能 。 

4.6 GATE. GATE 针对 生存 数据 利用 修正 的 泊 松 对 数 线性 混合 效应 模型 构建 似 然 函数 ， 
以 应 对 严重 的 删 失 率 ， 并 对 得 分 统计 量 进行 鞍点 近似 外 。 该 方法 也 能 够 校正 种 群 结构 和 亲 
缘 关 系 。 整 体 框架 与 SAIGE 类 似 。 

4.7 REGENIE. REGENIE 能 够 分 析 二 类 分 结局 的 极端 不 平衡 数据 ， 校 正方 法 包括 : 近 


似 Firth 校正 和 鞍点 近似 BJ。 相 较 于 SAIGE，REGENIE 不 再 估计 GRM， 而 是 采用 机 器 学 习 
的 方式 构造 种 群 结构 方差 ， 从 而 提高 运算 效率 。 为 了 得 到 快速 的 Firth 校正 结果 ， 该 方法 将 
协 变 量 效 应 6 作为 固定 值 进行 统计 量 构造 ， 从 而 避免 建立 联合 分 布 ， 以 提高 运算 效率 。 


五 、 讨 论 

大 型 生物 样本 库 和 全 基因 组 测序 数据 的 共享 促进 了 GWAS 数据 的 深度 挖掘 ， 也 推动 了 
遗传 统计 方法 的 发 展 。 本 文 主要 介绍 了 当前 GWAS 中 最 常用 的 两 种 处 理 极端 不 平衡 数据 的 
校正 方法 : Firth 校正 方法 和 鞍点 近似 方法 ， 并 总 结 了 相关 统计 软件 ， 骨 在 为 合理 分 析 


GWAS 数据 提供 方法 和 应 用 参考 。 

从 计算 速度 来 看 ， 为 实现 大 型 生物 样本 库 级 别 的 快速 分 析 ， 多 数 统计 方法 采用 得 分 检 
验 进行 统计 推断 ， 并 采用 校正 P 值 的 方法 实现 第 一 类 错误 的 控制 。 从 统计 性 能 来 看 ， 通 点 
近似 方法 能 够 较 好 地 控制 一 类 错误 ， 而 Firth 校正 方法 不 尽 如 和 人意 。 相 比 而 言 ， 参 数 估计 仍 
然 处 于 GWAS 分 析 的 核心 地 位 。 不 过 ， 研 究 者 需要 注意 以 下 三 点 : (0 SAIGE 等 采用 鞍点 
近似 的 方法 ， 只 调整 了 假设 检验 的 P 值 ， 未 校正 参数 估计 值 。 因 此 ，SNP 的 效应 和 标准 误 
的 估计 值 可 能 是 有 偏 的 。 不 建议 将 上 述 软件 所 得 结果 用 于 Meta 分 析 B1。(2) 无 论 是 Firth 校 
正 还 是 鞍点 近似 校正 的 统计 量 与 经 典 统计 量 完全 不 同 ， 包 括 : 效应 值 、 标 准 误 、 检 验 统计 
量 和 PP 值 。 因 此 ， 使 用 基于 汇总 数据 (summary-level data) 的 统计 方法 时 ， 需 要 事先 确认 可 行 
性 ， 包 括 但 不 限于 : 计算 SNP 遗传 度 、 构 建 多 基因 评分 (polygenetic risk score, PRS) 等 。(3) 
本 文中 所 列举 的 方法 是 均 针对 单位 点 (single-variant) 的 假设 检验 ， 获 得 的 系数 是 该 位 点 的 边 
际 效应 ， 因 此 不 宜 直接 将 回归 系数 用 于 构建 预测 模型 。 唯 一 可 用 于 构建 预测 模型 的 R 软件 
包 logistf 也 需要 在 Firth 校正 系数 的 基础 上 额外 校正 截 距 项 Bo。 
目前 极端 不 平衡 数据 的 GWAS 统计 分 析 已 日 臻 完善。 除了 前 文 所 提 及 用 于 二 分 类 和 生 
存 数据 的 方法 ， 鞍 点 近似 方法 还 可 拓展 至 多 分 类 结局 分 析 (POLMMBG”)、 基 因 - 环 境 交 互 作 


im 


H 


SAIGE-GENEP6., SAIGE-GENE4P7)), Æ 


| 4} Wr(SPAGEP?]) EJ 3 [A] f& (gene-based)4) Pr (fastSPA-GENE“4], POLMM-GENEP?!, 
g 


m} 
X 


音 的 统计 方法 满足 了 各 类 GWAS 数据 分 


析 的 需求 ， 实 现 了 更 合理 的 数据 分 析 ， 为 正确 鉴定 新 的 生物 标志 物 提供 了 理论 保障 和 
技术 支持 。 


参考 文献 


[1] 


[12] 


ZHOU W, NIELSEN J B, FRITSCHE L G, et al. Efficiently controlling for case-control 
imbalance and sample relatedness in large-scale genetic association studies [J]. Nat Genet, 
2018, 50(9): 1335-41. 


DENNY J C, BASTARACHE L, RITCHIE M D, et al. Systematic comparison of phenome- 
wide association study of electronic medical record data and genome-wide association 
study data [J]. Nat Biotechnol, 2013, 31(12): 1102-10. 


DEY R, ZHOU W, KIISKINEN T, et al. Efficient and accurate frailty model approach for 
genome-wide survival association analysis in large-scale biobanks [J]. Nat Commun, 2022, 
13(1): 5437. 


MA C, BLACKWELL T, BOEHNKE M, et al. Recommended joint and meta-analysis 
strategies for case-control association testing of single low-count variants [J]. Genet 
Epidemiol, 2013, 37(6): 539-50. 


DAI X, FU G, ZHAO S, et al. Statistical Learning Methods Applicable to Genome-Wide 
Association Studies on Unbalanced Case-Control Disease Data [J]. Genes (Basel), 2021, 
12(5). 


LEHMANN EL, ROMANO J P, CASELLA G. Testing statistical hypotheses [M]. Springer, 
1986. 


BERGER R L, CASELLA G. Statistical inference [M]. Duxbury, 2001. 


韩 栋 , 陈 征 , 陈 平 雁 , et al. 轮廓 似 然 函 数 及 其 应 用 [J]. 中 国 卫 生 统计 ,2012, 29(04): 478- 
80+83. 


EFRON B, HINKLEY D V. Assessing the accuracy of the maximum likelihood estimator: 
Observed versus expected Fisher information [J]. Biometrika, 1978, 65(3): 457-83. 


KENDALL M, STUART A, ORD J, et al. Vol. 2A: Classical inference and the linear model 
[J]. London [etc]: Arnold [etc], 1999. 


DEY R, SCHMIDT E M, ABECASIS G R, et al. A Fast and Accurate Algorithm to Test for 
Binary Phenotypes and Its Application to PheWAS [J]. Am J Hum Genet, 2017, 101(1): 37- 
49. 


DEY R, NIELSEN J B, FRITSCHE L G, et al. Robust meta-analysis of biobank-based 
genome-wide association studies with unbalanced binary phenotypes [J]. Genet Epidemiol, 
2019, 43(5): 462-76. 


[13] 
[14] 


[15] 


[16] 


[17] 


[18] 


[19] 


[20] 


[24] 


[25] 


[26] 


[27] 


[28] 


[29] 


COX D R. Partial likelihood [J]. Biometrika, 1975, 62(2): 269-76. 
MOORE D F. Applied survival analysis using R [M]. Springer, 2016. 


COX D R, SNELL E J. A general definition of residuals [J]. Journal of the Royal Statistical 
Society: Series B (Methodological), 1968, 30(2): 248-65. 


KANG H M, SUL J H, SERVICE S K, et al. Variance component model to account for 
sample structure in genome-wide association studies [J]. Nat Genet, 2010, 42(4): 348-54. 


YANG J, LEE S H, GODDARD M E, et al. GCTA: a tool for genome-wide complex trait 
analysis [J]. Am J Hum Genet, 2011, 88(1): 76-82. 


LOH P R, TUCKER G, BULIK-SULLIVAN B K, et al. Efficient Bayesian mixed-model 
analysis increases association power in large cohorts [J]. Nat Genet, 2015, 47(3): 284-90. 


BRESLOW N E, CLAYTON D G. Approximate inference in generalized linear mixed 
models [J]. Journal of the American statistical Association, 1993: 9-25. 


GILMOUR A R, THOMPSON R, CULLIS B R. Average information REML: an efficient 
algorithm for variance parameter estimation in linear mixed models [J]. Biometrics, 1995: 
1440-50. 


BI W, FRITSCHE L G, MUKHERJEE B, et al. A Fast and Accurate Method for Genome- 
Wide Time-to-Event Data Analysis and Its Application to UK Biobank [J]. Am J Hum 
Genet, 2020, 107(2): 222-33. 


FIRTH D. Bias reduction of maximum likelihood estimates [J]. Biometrika, 1993, 80(1): 
27-38. 


GREENLAND S, MANSOURNIA M A. Penalization, bias reduction, and default priors in 
logistic and related categorical and survival regressions [J]. Statistics in medicine, 2015, 
34(23): 3133-43. 


HEINZE G. The application of Firth’s procedure to Cox and logistic regression [M]. 
Technical report 10. Department of Medical Computer Sciences, Section of Clinical 
Biometrics .... 1999. 


HEINZE G, SCHEMPER M. A solution to the problem of separation in logistic regression 
[J]. Stat Med, 2002, 21(16): 2409-19. 


LUGANNANI R, RICE S. Saddle point approximation for the distribution of the sum of 
independent random variables [J]. Advances in applied probability, 1980, 12(2): 475-90. 


FEUERVERGER A. On the empirical saddlepoint approximation [J]. Biometrika, 1989, 
76(3): 457-64. 


HOSMER JR D W, LEMESHOW S, STURDIVANT R X. Applied logistic regression [M]. 
John Wiley & Sons, 2013. 


HOSMER JR D W, LEMESHOW S, MAY S. Applied survival analysis: regression 
modeling of time-to-event data [M]. John Wiley & Sons, 2008. 


[30] 


[31] 


[32] 


[33] 


[34] 


[35] 


[36] 


[37] 


PUHR R, HEINZE G, NOLD M, et al. Firth's logistic regression with rare events: accurate 
effect estimates and predictions? [J]. Stat Med, 2017, 36(14): 2302-17. 


MBATCHOU J, BARNARD L, BACKMAN J, et al. Computationally efficient whole- 
genome regression for quantitative and binary traits [J]. Nat Genet, 2021, 53(7): 1097-103. 


BI W, ZHOU W, DEY R, et al. Efficient mixed model approach for large-scale genome- 
wide association studies of ordinal categorical phenotypes [J]. The American Journal of 
Human Genetics, 2021, 108(5): 825-39. 


BI W, ZHAO Z, DEY R, et al. A fast and accurate method for genome-wide scale phenome- 
wide Gx E analysis and its application to UK biobank [J]. The American Journal of Human 
Genetics, 2019, 105(6): 1182-92. 


ZHAO Z, BI W, ZHOU W, et al. UK Biobank whole-exome sequence binary phenome 
analysis with robust region-based rare-variant test [J]. The American Journal of Human 
Genetics, 2020, 106(1): 3-12. 


BI W, ZHOU W, ZHANG P, et al. Scalable mixed model methods for set-based association 
studies on large-scale categorical data analysis and its application to exome-sequencing 
data in UK Biobank [J]. Am J Hum Genet, 2023, 110(5): 762-73. 


ZHOU W, ZHAO Z, NIELSEN J B, et al. Scalable generalized linear mixed model for 
region-based association tests in large biobanks and cohorts [J]. Nature genetics, 2020, 
52(6): 634-9. 


ZHOU W, BI W, ZHAO Z, et al. SAIGE-GENE- improves the efficiency and accuracy of 
set-based rare variant association tests [J]. Nat Genet, 2022, 54(10): 1466-9. 


